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This program is an important tool in the study of alternative routes between the Earth and 
the Moon. Dr. John Aired was the NASA Technical Monitor for the contract under which this 
program was produced. Mr. Andy Petro was the NASA Task Manager for this particular task. 
Mr. Bill Stump was the Eagle Project Manager for the contract under which this program was 
produced. The program was written by Jack Funk, originally in Quick Basic, and translated into 
Fortran by Mr. Bill Engblom. Mr. Engblom also prepared the documentation. 
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1.0 Introduction 


The program LIB RATE calculates velocities for trajectories from low earth orbit (LEO) 
to four of the five libration points (L2, L3, L4, L5), and from low lunar orbit (LLO) to libration 
points LI and L2. Librations points (LP) are defined as locations in space that orbit the Earth 
such that they are always stationary with respect to the Earth-Moon line. Libration point #2 (L2) 
is located between the Earth and Moon where the gravitational attraction from both bodies are 
equal. LI and L3 are located behind the Moon and Earth, respectively, such that the pull of the 
Earth and Moon together just cancel the centrifugal acceleration associated with the libration 
point’s orbit. L4 and L5 are located half-way between the Earth and the Moon and 60° off the 
Earth-Moon line to the left and right, respectively. Hence, the Earth, Moon, and all libration 
points, lie in the same plane. 

The flight to be analyzed departs from a circular orbit of any altitude and inclination 
about the Earth or Moon and finishes in a circular orbit about the Earth at the desired libration 
point within a specified flight time. First, the departure orbit is made into a more eccentric orbit 
(ellipse or hyperbola) with an initial AV in order to reach the libration point while meeting the 
flight time constraint. The less the desired flight time, the more eccentric the orbit, and the 
larger the initial AV required. The least eccentric elliptic orbit would require the minimum AV 
and the maximum flight time. A second AV is then needed once the elliptic or hyperbolic flight 
path has reached the libration point in order to change the velocity vector of the eccentric 
trajectory to that of the libration point’s orbit (circularize). So, the more eccentric the orbit, the 
larger the velocity change. This second bum must also account for the inclination of the 
eccentric trajectory with respect to the Earth-Moon-LP plane. 
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This program produces a matrix of the AV’s needed to complete the desired flight. The 
user specifies the departure orbit (location and altitude), and the maximum flight time. A matrix 
is then developed with 10 inclinations (with respect to the Earth-Moon-LP plane), ranging from 
0° to 90°, forming the columns, and 19 possible flight times, ranging from the flight time (input) 
to 36 hours less than the input value, in decrements of 2 hours, forming the rows. This matrix is 
presented in three different reports including the total AV’s, and both of the AV components 
discussed above. 

Section 2.0 of this document describes the input required from the user to define the 
flight. Section 3.0 describes the contents of the three reports that are produced as outputs. 
Section 4.0 includes the instructions needed to execute the program. 

A more detailed description of the process used in LEBRATE has been included as 
Appendix D (main program), Appendix E (in-program subroutine), and Appendices F, G, H, and 
I (external subroutines). LEBRATE was derived in part from the PLANECHG program (also 
produced under this contract), discussed in a different, more detailed documentation report. 
Therefore, for a more in-depth look at many of the equations, variables, and conventions used in 
LIBRATE, please consult the PLANECHG documentation. 
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2.0 Program Inputs 


The following paragraphs discuss the inputs provided by the user. The prompt is the 
message displayed by the program onto the screen. The input variable is the variable assigned to 
the user’s response. The description provides information about how to respond to the prompt. 

1 . Prompt: INPUT EARTH OR MOON 

Input variable:BODY 

Description: Enter either MOON or EARTH as a departure orbit location. EARTH is 

the default value. 

2. Prompt: INPUT PERIGEE ALTITUDE (NMI) 

Input variable:HPE 

Description: Enter departure orbit altitude above the surface of the departure location 

(Earth or Moon), in nautical miles. A circular orbit is assumed. 

3. Prompt: INPUT LIB RATION POINT NUMBER 

Input variable:NLP 

Description: Enter the number (1, 2, 3, 4 or 5) of the desired libration point. Remem- 

ber that LI (behind the Moon), L2 (between Earth and Moon), and L3 
(behind Earth) are all on the Earth-Moon line. L4 and L5 are located 
midway between the Earth and Moon and 60° off the Earth-Moon line to 
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the left and right, respectively. Recall that LEBRATE cannot calculate 
trajectories from the Earth to LI or the Moon to L3, L4, or L5. Other 
programs produced under this contract do these calculations, however. 

3. Prompt: INPUT FLIGHT TIME 

Input variable:FLTIM 

Description: Enter the flight time allowed for the transfer, in hours. If this value is 

larger than the calculated maximum flight time for such a trajectory then a 
message will appear indicating the maximum flight time allowed followed 
by another flight time input prompt. This may occur several times 
because the program approximates the maximum flight time based on the 
flight time input. Simply continue to enter flight times less than those 
time constraints issued by the program until a value is accepted (no error 
message). 
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3.0 Program Outputs 


This section describes the contents of each of the three reports generated by the program. 
These reports may be found in the output file, LBBRATE.OUT. Samples of report #1, #2, and #3 
have been included herein, starting on the following page. 

Report #1: Total Delta Velocity for Earth (or Moon) Transfer Trajectory to Libration Point 

This report is a matrix of total delta velocities in ft/sec that are needed to complete the 
transfer from the departure location (Earth or Moon) to the desired libration point. Each cell 
corresponds to a particular flight time and inclination (of departure orbit with respect to Earth- 
Moon-LP plane). 

Report #2: Delta Velocity at Libration Point for Earth (or Moon) Transfer Trajectory to 
Libration Point 

This report is a matrix of delta velocities in ft/sec that are needed at the libration points to 
correct the velocity vectors of the eccentric trajectories to that of the libration point’s orbit. Each 
cell corresponds to a particular flight time and inclination (of departure orbit with respect to 
Earth-Moon-LP plane). 

Report #3: Delta Velocity at Earth (or Moon) Orbit for Transfer Trajectory to Libration Point 

This report is a matrix of delta velocities in ft/sec that are needed to make the circular 
departure orbit into a trajectory (ellipse or hyperbola) that will reach the libration point while 
meeting the flight time constraint. Each cell corresponds to a particular flight time and inclina- 
tion (of departure orbit with respect to Earth-Moon-LP plane). 



OF TOTAL VELOCITY NEEDED FOR EARTH TRANSFER TRAJECTORY TO LIBRATION POINT #2 

RUN DATE 23-SEP-88 RUN TIME 13:47:50 
CIRCULAR ORBIT ALTITUDE 250.0 
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OF DELVEL AT LIBRATION POINT FOR EARTH TRANSFER TRAJECTORY TO LIBRATION POINT #2 

RUN DATE 23-SEP-88 RUN TIME 13:47:50 
CIRCULAR ORBIT ALTITUDE 250.0 
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4.0 Program Execution Instructions 


The following instructions describe the steps to be taken by the user to execute this program: 

A. Obtain access to the DEC VAX minicomputer and sign on with user identification. 

B. At the $ prompt, type RUN LIB RATE . 

C. When prompted by the program, enter the program inputs. See section 2.0 for a 
discussion of the inputs. 

D. After the last input has been entered, the program will execute for approximately 1 
minute. Upon completion, the message FORTRAN STOP will appear, followed by 
the $ prompt. 

E. The program outputs will be placed in a file named U® RATE. OUT;### where ### is 
a system generated number of the report. To print the most recently generated report, 
type the following at the $ prompt: TYPE LIBRATE.OUT 

F. To re-execute the program with new parameters, begin again at step (B) above. 
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Appendix A. Program Flow Chart 
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Appendix B. Code Listing 


oo oooooo 


0 ************************************************************* 
C *** Libration Point (LP) Program *** * 

C ************************************************************* 

C * from Earth to Points 2, 3, 4, and 5 and * 

C * from Moon to Points 1 and 2 * 

q ************************************************************* 
C * Written in Quick Basic by: Jack Funk * 

C * Translated into Fortran and documented by: Bill Engblom * 

0 ************************************************************* 
C 

C MAIN 

C 

IMPLICIT REAL*16 (A-H,0-Z) 

CHARACTER* 5 TRAJE, BODY 
CHARACTER* 3 2 HEAD1 , HEAD 2 , HEAD 3 
CHARACTER* 10 TIMP, DATP 


DIMENSION XXL (5) , Y YL ( 5 ) , AINCEP(IO), FLTIMP(19) 
DIMENSION PAGE1 (15, 19) , PAGE2 (15,19), PAGE3(15,19) 

OPEN OUTPUT FILE 

OPEN (UNIT = I, FILE = ' LIBRATE . OUT ' f STATUS = ' NEW' ) 

OUTPUT TO FILE; IP: OUTPUT TO SCREEN; IS 

IP =1 

IS =5 

II =1 

JJ =1 

NS TOP = 1 

DPR = 57.29578 

PI * 3.141593 

CMUE = 1 . 407647E+16 

CMUM = 1 . 731400E+14 

FTNM = 6076.115 

REE = 20925741. 

REMO = 5 . 7039E+6 
FTM = .3048 

OUTBOUND TRAJECTORIES (MODE = 1) 

MODE HAS BEEN CHANGED TO MD 
MD - 1 

CALL DATE (DATP) 

CALL TIME (TIMP) 

PRINT *,' INPUT EARTH OR MOON ' 

READ (6,5) BODY 
5 FORMAT (A5) 

IF (BODY .NE. 'MOON') BODY = 'EARTH' 

PRINT *, 'INPUT PERIGEE ALTITUDE (NMI) ' 

READ *, HPE 



RPE = HPE * FTNM + REE 

PRINT *, ' INPUT LIBRATION POINT NUMBER ' 

READ * , NLP 
AINCEO = 0. 

AINCE = AINCEO / DPR 
20 CONTINUE 

PRINT * , ' INPUT FLIGHT TIME (HR) ' 

READ *, FLTIM 
IF (FLTIM .EQ. 0) GOTO 20 
FTIM = FLTIM 
25 CONTINUE 

C CALCULATE EARTH-MOON DISTANCE 

CALL MOON (T, RAM, DECM, RM) 

CALL POSVELMO (TIEM, RREMER, XDLO, YDLO) 

RREM - RM * REE 
C WRITE (IS, 30) 

C 30 FORMAT (' FTIM DVTOT DVCIR DELV VELX VXX VYX 
C * VXLP VYLP RREM' ) 

C CALCULATE LIBRATION POINT LOCATIONS, MAX FLIGHT TIME 

I CALL = 1 
GOTO 70 
1000 CONTINUE 
C 

35 CONTINUE 
C 

C UPDATE LIBRATION POINT LOCATIONS 

I CALL = 2 
GOTO 70 
2000 CONTINUE 
40 CONTINUE 

C TIMEE HAS BEEN CHANGED TO TIEM 

TIEMS = TIEM 

CALL GAMACALC (RPE , VELE, RRL, CMUE, CSGAME, VPE, VCIRE, 

* TIEM1, TRAJE, THETAE) 

VELE2 = VELE + 1. 

CALL GAMACALC (RPE, VELE2, RRL, CMUE, CSGAME, VPE, VCIRE, 

* TIEM2, TRAJE, THETAE) 

DELVELE = 1 . / (TIEM2 - TIEM1) * (FTIM - TIEM1 - TIEMM) 

IF (DELVELE .LT. 500.) THEN 
VELE = VELE + DELVELE 
ELSE 

VELE = VELE + DELVELE/QABS (DELVELE) * 500. 

END IF 

IF (VELE .LT. VPMIN) VELE = VPMIN 

CALL GAMACALC (RPE, VELE, RRL, CMUE, CSGAME, VPE, VCIRE, 

* TIEM, TRAJE, THETAE) 

C 

IF (TIEM .EQ. 0.) GOTO 60 
TIEMT = TIEM 


VZX 



IF (QABS (FTIM - TIEMT) .GT. 1.) GOTO 35 
C 

C MODE HAS BEEN CHANGED TO MD 

GAMAE = FLOAT (MD) * QATAN (QSQRT (1 . - CSGAME**2 . ) /CSGAME) 

ALONE = 180. /DPR + ALONX 

ALATE = QATAN (ZZLP/QSQRT (XXLP**2 + YYLP**2) ) 

50 CONTINUE 

C CALCULATION OF EARTH ORBIT AZIMUTH AT SPHERE OF INFLUENCE 

IF (AINCE .NE. 0.) THEN 

PHIE = PI-QASIN (ZZ/ (RRL*QSIN (AINCE) ) ) 

SAZM = QCOS (AINCE) / QCOS (ALATE) 

CAZM = QSIN (AINCE) * QCOS (PHIE) / QCOS (ALATE) 

AZME = QATAN 2 (SAZM, CAZM) 

ELSE 

AZME = PI/2. 

END IF 
C 

CALL VELTRANS (VELE , GAMAE , AZME , ALATE , ALONE , VXE , VYE , VZE , VELE ) 
DELVEL = QSQRT ( (VXE+OMEGM*YYLP) **2 . + (VYE-OMEGM* (XXLP 
* - XXBC) ) **2 . + VZE**2 . ) 

DVCIRE = VPE - VCIRE 
DVTSAV = DVTOTAL 
DVTOTAL = DELVEL + DVCIRE 
C 

C CALCULATE POSITION AND VELOCITY OF MOON 

C ITTERATE FOR FLIGHT TIME TO SPHERE OF INFLUENCE 

C 

TIM = (TIMJ-2451545 . ) /36525 . +TIEM/ 876600 . 

CALL POSVELMO (TIM, RREMER, XDLO, YDLO) 

RREM = RREMER* REE 

IF (QABS (TIEM-TIEMS) .GT. .09) GOTO 35 
60 CONTINUE 

C WRITE (IS, 65) TIEM, DVTOTAL, DVCIRE, DELVEL, VELE, VXE, VYE, VZE, 

C * -OMEGM*YYLP, OMEGM* (XXLP-XXBC) , RREMER 

C 65 FORMAT (F5 . 1, 2X, F6.0,1X,F6.0, IX, F5 . 0, 2X, F5 . 0 , IX, 3 (F6 . 0 , IX) 

, IX, F6 . 0, 

C * IX, F6 . 0, 2X, F6 . 3) 

AINCEP (II) = AINCE * DPR 
FLTIMP(JJ) = FTIM 
PAGE1 (II, J J) = DVTOTAL 
PAGE2 (II, J J) = DELVEL 
PAGE3 (II, J J) * DVCIRE 
FTIM = FTIM - 2. 

JJ = JJ + 1 

NSTOP = NSTOP + 1 

IF (NSTOP .LT. 20) GOTO 40 

FTIM = FLTIM 

JJ = 1 

NSTOP = 1 



AINCE = AINCE + 10. /DPR 
II = II + 1 

IF (AINCE .LT. 91. /DPR) GOTO 40 
HEAD1 = ' MAP OF TOTAL VELOCITY NEEDED' 

CALL LPPAGE (PAGE1, BODY, HEAD1, NLP, DATP, TIMP, HPE, AINCEP, 

* FLTIMP) 

HEAD 2 = 'MAP OF DELVEL AT LIBRATION POINT' 

CALL LPPAGE (PAGE2, BODY, HEAD2, NLP, DATP, TIMP, HPE, AINCEP, 

* FLTIMP) 

IF (BODY .EQ. 'EARTH') THEN 

HEAD 3 = 'MAP OF DELVEL AT EARTH ORBIT' 

ELSE 

HEAD 3 = 'MAP OF DELVEL AT LUNAR ORBIT' 

END IF 

CALL LPPAGE (PAGE3, BODY, HEAD3, NLP, DATP, TIMP, HPE, AINCEP, 

* FLTIMP) 

GOTO 3000 

70 CONTINUE 

C CALCULATION OF LIBRATION POINT POSITIONS 

RREM = RREMER * REE 
XXBC = .01215052 * RREM 
OMEGM = YDLO/ (RREM - XXBC) 

XXL (1) = 1.15567 * RREM 
YYL(l) = 0. 

XXL (2) = 0.83691 * RREM 
YYL ( 2 ) = 0. 

XXL (3) = -1.005062 * RREM 
YYL (3) = 0. 

XXL (4) = 0.5 * RREM 

YYL (4) - -RREM * QSIN (60 . /DPR) 

XXL (5) = 0.5 * RREM 

YYL (5) = RREM * QSIN ( 60 . /DPR) 

ZZLP = 0. 

C CHANGE DATA FROM EARTH TO MOON 

IF (BODY .EQ. 'MOON') THEN 
RPE = HPE*FTNM + REMO 
XXLP = XXL (NLP) - RREM 
CMUE - CMUM 
ELSE 

XXLP = XXL (NLP) 

END IF 

C CALCULATE Y POSITION, RADIUS, AND LONGITUDE OF LIBRATION POINT 

NLP 

YYLP = YYL (NLP) 

RRL = QSQRT (XXLP* *2 . + YYLP**2. + ZZLP**2.) 

ALONX = QATAN2 (YYLP, XXLP) 

C CALCULATION OF MIN VEL AT LIB POINT FOR ORBIT CONTAINING 

C RRL AND RPE 

AAE = (RRL + RPE) /2 . 



VPMIN = QSQRT (CMUE/AAE*RPE/RRL) * 1.001 
IF (VELE .LT. VPMIN) VELE = VPMIN 

CALL GAMACALC (RPE , VPMIN, RRL, CMUE , CSGAME , VPE , VCIRE , AMAXFT , 

* TRAJE, THETAE) 

IF (AMAXFT .LT. FTIM) THEN 

PRINT FLIGHT TIME IS GREATER THAN MAX FLIGHT TIME ' , 

* AMAXFT 
GOTO 20 

END IF 

IF (ICALL. EQ. 1) GOTO 1000 
IF (ICALL .EQ. 2) GOTO 2000 
3000 STOP 
END 


SUBROUTINE LPPAGE (PAGE , BODYP , HEADP , NLP , DATP , TIMP , HPE , AINCEP , 

* FLTIMP) 

IMPLICIT REAL*16 (A-H,0-Z) 

CHARACTER* 5 BODYP 
CHARACTER* 10 TIMP, DATP 
CHARACTER* 32 HEADP 

DIMENSION PAGE (15,19) , AINCEP(IO), FLTIMP(19) 

IP = 1 
IS = 5 

WRITE (IP, 5) CHAR (12) 

5 FORMAT (' ' , Al) 

WRITE (IP, 10) HEADP, BODYP, NLP 
10 FORMAT (T2, A32, ' FOR ',A5,' TRANSFER TRAJECTORY 

* TO LIBRATION POINT #',I1) 

WRITE (IP, 20) DATP, TIMP 

20 FORMAT (T25,'RUN DATE ' ,A9, ' RUN TIME ',A8) 

WRITE (IP, 25) HPE 

25 FORMAT (T25, ' CIRCULAR ORBIT ALTITUDE ',F6.1) 

WRITE (IP, 30) AINCEP (1), AINCEP (2), AINCEP (3), AINCEP (4) , 

* AINCEP (5) , AINCEP (6) , AINCEP (7) , AINCEP (8) , AINCEP (9) , AINCEP (10) 

30 FORMAT (/, ' INCL>' , 10 (2X, F5 . 1) , /, ' FLIGHT',/,' TIME (HR)') 

DO 40 NPI = 1,19 

WRITE (IP, 35) FLTIMP (NPI) , PAGE (1, NPI) , PAGE (2, NPI) , 

* PAGE (3, NPI) , PAGE (4, NPI) , PAGE (5, NPI) , PAGE (6, NPI) , PAGE (7, NPI) , 

* PAGE (8, NPI) , PAGE (9, NPI) , PAGE (10, NPI) 

35 FORMAT (IX, F6 . 1 , 10 ( IX, F6 . 0) ) 

40 CONTINUE 
RETURN 
END 



SUBROUTINE POSVELMO (TIM, RRM, XDLO, YDLO) 
IMPLICIT REAL * 16 <A-H,0-Z) 

DELT = 0. 5/36525. /24./3600. 

T1 = TIM-DELT 

CALL MOON (Tl, RAM, DECM,RMl) 

T2 = TIM+DELT 

CALL MOON (T2, RAM, DECM,RM2) 

XDLO = (RM2-RM1) *20925741 . 

RRM = (RM2+RM1) /2 . 

RRM = RRM 

RRMB = RRM -7 . 41278 9E-01 
YDLO = 200570 .2/RRMB 
RETURN 
END 


SUBROUTINE MOON (T, RAM, DECM, RM) 

C FINDS LOCATION OF MOON IN EQUATORIAL COORDS. AT ANY TIM 

C REF: '87 ASTRONOMICAL ALMANAC 

C T IS JULIAN CENTURIES SINCE YEAR 2000 

C LAM IS MOON'S ECLIPTIC LONGITUDE 

C BETA IS MOON'S ECLIPTIC LATITUDE 

C PIE IS HORIZONTAL PARALLAX 

C RM IS DIST. TO MOON IN EARTH RADII 

C RAM IS RT. ASCENSION OF MOON 

C DECM IS MOON'S DECLINATION 

C SD IS SEMIDIAMETER OF MOON'S ORBIT 

IMPLICIT REAL * 16 (A-Z) 

C PRINT *, ' MOON' 

P = 3.1415926535 

C = P / 180. 

LAM = C*218 . 32+C*481267 . 883*T+C* 6.29 * QSIN(C * 134.9 + C * 
*477198.85 * T) - C * 1.27 * QSIN(C * 259.2 - C * 413335.38 * 
*T) + c * .66 * QSIN (C * 235.7 + c * 890534. 23*T) 

LAM = LAM + C * .21 * QSIN(C * 269.9 + C * 954397. 7*T) - C * 
*.19 * QSIN (C * 357.5 + C * 35999.05 * T) - C * .11 * 

*SIN (C * 186.6 + C * 966404. 05*T) 

beta « c*5 . 13*QSIN (c*93 . 3 + c * 483202.03 * T) + c * 

* .28 * QSIN (C * 228.2 + C * 960400.87 * T) -c*.28*QSIN 
* (c*318 . 3+c* 60003. 18*T) -c* . 17*QSIN (c*217 . 6-c*407332 . 2 * T) 

pie = c* . 9508+c* . 0518*COS (c*134 . 9 + c * 477198.85 * T) + 



*c*.0095*QCOS(c*259.2-c*413335.38*T)+c* .0078*COS(c * 23 
*5 . 7+c*890534 . 23*T) +c* . 0028*QCOS (c*269 . 9+c*954397 . 7* T) 

SD = .2725 * pie 
RM = 1. / QSIN(pie) 

1 = QCOS (beta) * QCOS (LAM) 

M = .9175 * QCOS (beta) * QSIN(LAM) - .3978 * QSIN(beta) 

n = .3978 * QCOS (beta) * QSIN(LAM) + .9175 * QSIN(beta) 

RAM = QATAN2 (M, 1) 

DECM = QASIN(n) 

RETURN 

END 


SUBROUTINE GAMACALC (RPX, W, RRX, CMUX, COSGAMX, VPX, VCIRX, 
*TIMX, TRAJ, THE TAX, DPR, FTNM) 

IMPLICIT REAL * 16 (A-H, 0-Z) 

CHARACTER* 5 TRAJ 
IHYPER = 1 
TRAJ = 'ELIPT' 

QQX = RRX* W** 2. /CMUX 

IF (QQX-2 .LT. 1.0E-06) QQX = QQX - 1.0E-06 
IF (QQX .GT. 2.) THEN 
IHYPER = -1 
TRAJ = ' HYPER' 

C PRINT *,' ***** TRAJECTORY IS HYPERBOLIC ' 

END IF 

AAX = RRX/ (2. -QQX) 

IF (AAX. GT . 1.0E12 .OR. AAX. LT. -1.0E12) AAX = -1.0E12 

EEX = 1 . — RPX/AAX 

PPX = AAX* (1 . -EEX ** 2.) 

COSGAMX = QSQRT (RPX/RRX* (1 . +EEX) /QQX) 

GAMAX = QACOS (COSGAMX) 

VPX = QSQRT (CMUX* (1 . +EEX) /RPX) 

VCIRX = QSQRT (CMUX/RPX) 

COSTHETAX = (PPX/RRX-1 . ) /EEX 

THETAX = QACOS (COSTHETAX) 

COSAEX = (EEX+COSTHETAX) / (1 . +EEX*COSTHETAX) 

ERRRP ■ RRX -AAX* (1+EEX) 

IF (ERRRP .GT. 0. .AND. AAX . GT . 0.) THEN 

WRITE (IS, 407) ERRRP/ 6076.1155, ALAT*DPR, ALON*DPR 
407 FORMAT (' RADIUS > APOGEE BY NMI ',F7.5, F8.0, F7.5) 
WRITE (IS, 417) QQX, AAX/FNTM, EEX, GAMAX*DPR, VPX 
417 FORMAT (' QQX = ',F7.5,' AAX = ',F8.0,' EEX = ',F7.5, 

* ' GAMAX = ' , F5.1,' VPX = ' ,F7 .1) 

ENDIF 



c 


IF(QQX .GT. 2.) GOTO 101 

CALC FLIGHT TIME FOR ELIPTICAL ORBITS 
IF (ERRRP .GT. 0.) THEN 
TIMX = 0. 

GOTO 199 
END IF 

AEX = QACOS (COSAEX) 

SINAEX = QSQRT(1 -COSAEX* *2.) 

TIMX * QSQRT (AAX ** 3. / CMUX) * (AEX-EEX* SINAEX) /3600 . 
GOTO 199 
101 CONTINUE 

C CALC FLIGHT TIME FOR HYPERBOLIC ORBITS 

CO SHF * COSAEX 
SINHF = QSQRT (COSHF**2 . -1 . ) 

FFX = QLOG (COSHF+QSQRT (COSHF**2 . -1 . ) ) 

TIMX = QSQRT (-1 . *AAX*AAX*AAX/CMUX) * (EEX* SINHF -FFX) /3600 . 
199 CONTINUE 
RETURN 
END 


SUBROUTINE VELTRANS (VEL, GAMA, AZM, ALAT, ALON, VXX, VYX, VZX, VMAG) 
IMPLICIT REAL * 16 (A-Z) 

RRD = VEL * QSIN (GAMA ) 

RPHID = VEL * QCOS (GAMA ) 

VLON « RPHID * QSIN (AZM ) 

VLAT = RPHID * QCOS (AZM ) 


VXR 

= -RRD 

* 

QCOS 

(ALAT 

) 

* 

QCOS 

(ALON 

) 

VYR 

= -RRD 

* 

QCOS 

(ALAT 

) 

* 

QSIN 

(ALON 

) 

VZR 

= RRD 

* 

QSIN 

(ALAT 

) 





VXLA 

= VLAT 

* 

QSIN 

(ALAT 

) 

* 

QCOS 

(ALON 

) 

VYLA 

- VLAT 

* 

QSIN 

(ALAT 

) 

* 

QSIN 

(ALON 

) 

VZLA 

= VLAT 

* 

QCOS 

(ALAT 

) 





VXLO 

= VLON 

* 

QSIN 

(ALON 

) 





VYLO 

= -VLON 

if 

QCOS 

(ALON 

) 





VZLO 

II 

o 

• 

o 









VXX 

= VXR + 

VXLA + 

VXLO 






VYX 

= VYR + 

VYLA + 

VYLO 






VZX 

= VZR + 

VZLA + 

VZLO 






VMAG 

= QSQRT 

(VXX **2 . + VYX 

** 2.+ 

VZX ** 


RETURN 

END 



Appendix C. Program Variables 


INPUT 

VARIABLE DESCRIPTION 

AINCEO Inclination of departure orbit (LEO or LLO) with respect to the Earth-Moon-LP 
plane (degrees) 


BODY Body of departure (’EARTH’ or ’MOON’) 

FLTIM Flight time constraint for trajectory (hours) 

HPE Holding perigee altitude of departure orbit (nautical miles) 

CONSTANT VALUE DESCRIPTION 

C 

11/180 

Degrees per radian (deg./rad.) 

CMUE 

1.407647E+16 

Gravitational parameter of the Earth (ftVsec 2 ) 

CMUM 

1.731432E+14 

Gravitational parameter of the Moon (ft 3 /sec 2 ) 

DPR 

57.29578 

Degrees per radian (deg./rad.) 

FTM 

0.3048 

Meters per foot (m/ft) 

FTNM 

6,076.115 

Feet per nautical mile (ft/nm) 

MD 

’OUTBOUND’ 

Leg of trip on which to perform calculations (outbound) 

PI 

3.141593 

11 (dimensionless) 

REE 

20,925,741 

Radius of the Earth (ft) 

REMO 

5,703,900 

Radius of the Moon (ft) 


VARIABLE 

AAE 

AAX 

AEX 

AINCE 

AINCEP 

ALAT 

ALATE 

ALATMIN 

ALATSIM 

ALON 

ALONE 

ALONMIN 

ALONSIM 

ALONX 

AMAXFT 

AZM 

AZME 

BETA 

CAZM 

CMUX 


DESCRIPTION 

Semi-major axis of least eccentric Earth (or Moon)-to-LP trajectory. Perigee at 
departure orbit, apogee at libration point 

Semi-major axis of one of the transfer orbits (Gamacalc Subroutine) 

Eccentric anomaly of one of the transfer orbits (Gamacalc Subroutine) 

Departure orbit inclination (rad.) with respect to Earth-Moon-LP plane 
Array of departure orbit inclinations (Lppage Subroutine) 

Latitude of LP point, measured from Earth-Moon plane, always 0° (Gamacalc 
Subroutine) 

Latitude of LP, measured from Earth-Moon plane, always 0° 

Latitude of LP associated with the minimum total AV 

Latitude of LP associated with the minimum heading correction AV 

Longitude of LP, measured from Earth-Moon line zero longitude at departure 
body, starts from -X direction (Veltrans Subroutine) 

Longitude of LP, measured from zero longitude at Earth, starts from -X direction 

Longitude of LP associated with the minimum total AV 

Longitude of LP associated with the minimum heading correction AV 

Longitude of the LP, measured from the Earth-Moon line at the departure body, 
starts from -X direction 

Maximum flight time to LP from perigee altitude (least eccentric trajectory) 

Azimuth of the transfer orbit upon reaching the LP, measured clockwise from +Y 
direction of departure body (Veltrans Subroutine) 

Azimuth of the transfer orbit upon reaching the LP, measured clockwise from +Y 
direction of departure body 

Moon’s ecliptic latitude 

Cosine of the azimuth 

Earth or Moon gravitational parameter (Gamacalc Subroutine) 



COSAEX 


Cosine of the eccentric anomaly of one of the transfer orbits 

CSGAME Cosine of the flight path angle at LP of the Earth (or Moon)-to-LP trajectory 

COSGAMX Cosine of the flight path angle at LP of one of the transfer orbits (Gamacalc 
Subroutine) 

COSHF Hyperbolic cosine of the eccentric anomaly of one of the transfer orbits 
COSTHETAXCosine of the true anomaly of one of the transfer orbits 


DATP 

DECM 

DELT 

DELVEL 

DELVELE 

DVCIRE 

DVMIN 

DVSIM 

DVTOTAL 

DVTSAV 

ERRRP 

EEX 

FFX 

FLTIMP 

GAMA 

GAMAE 

GAMAX 

HEADP 


Today’s date 

Declination of the Moon (not used) 

A fraction of TEM that represents a half-second 
AV needed at the LP point to correct the velocity vector 

Extrapolated AV increment which is added to existing guess for the velocity at the 
LP of the Earth-to-LP trajectory (VELE) 

AV between departure circular orbit and Earth (or Moon)-to-LP trajectory 

Hold variable for minimum total AV 

Hold variable for minimum heading correction AV 

Total AV for flight 

Temporary storage for total AV 

Difference between orbital range at LP and apogee range. Orbital radius must be 
less than apogee radius or error message results 

Eccentricity of one of the transfer orbits 

Eccentric anomaly of one of the transfer orbits 

Array of flight times (Lppage Subroutine) 

Flight path angle at LP of Earth (or Moon)-to LP trajectory (Veltrans Subroutine) 
Flight path angle at LP of Earth (or Moon)-to LP trajectory 
Flight path angle of one of the transfer orbits (Gamacalc Subrou 
Heading for report (Lppage Subroutine) 



HEAD1 

HEAD2 

HEAD3 

IHYPER 

ICALL 

n 

IP 

IS 

JJ 

L 

LAM 

M 

N 

NSTOP 

OMEGM 

OMEGE 

PAGE 

PAGE1 

PAGE2 

PAGE3 

PHIE 

PIE 

PPX 


Heading for report #1 
Heading for report #2 
Heading for report #3 

Indicator that describes whether an orbit is hyperbolic 

Flag signifying where the program will resume after the internal 
subroutine (GOTO 70) has been completed 

Counter for rows of PAGE arrays 

Unit number for formatted writes to output file 

Unit number for formatted writes to the screen 

Counter for columns of PAGE arrays 

A geocentric direction cosine 

Moon’s ecliptic longitude 

A geocentric direction cosine 

A geocentric direction cosine 

Counter used to determine when to get a new inclination and reset the flight time 

Angular velocity of the Moon 

Argument of perigee of departure orbit 

Array of output matrix data (Lppage Subroutine) 

Array holding matrix of total AV data 

Array holding matrix of heading correction AV data 

Array holding matrix of departure AV data 

Angle between the Earth-to-Moon plane and the departure body-to-LP line 
(always 180°). Used to determine azimuth angle 

Horizontal parallax 

Semi-latus rectum of one of the transfer orbits 



QQX 

RAM 

RM 

RM1 

RM2 

RPE 

RPHID 

RPX 

RRD 

RREM 

RREMER 

RRL 

RRM 

RRMB 

RRX 

SAZM 

SD 

SINAEX 

SINHF 

T 

T1 

T2 

THETAE 

THETAX 


Vis-viva parameter for an orbit 
Right ascension of the Moon (not used) 

Distance from the Earth to the Moon in Earth radii at time T 
Distance from the Earth to the Moon in Earth radii at time T1 
Distance from the Earth to the Moon in Earth radii at time T2 
Distance from departure body center to perigee altitude 
Orbital path component of the velocity vector 

Distance from Earth or Moon center to perigee altitude (Gamacalc Subroutine) 

Radial component of the velocity vector 

Distance from the Earth to the Moon in feet 

Moon’s distance from the Earth in Earth radii 

Distance from the departure body center to the libration point 

Average distance from the Earth to the Moon, using RM1 and RM2 

Distance from the Earth-Moon baricenter to the Moon 

Distance from departure body (Earth or Moon) to LP 

Sine of the Azimuth 

Semi-diameter of the Moon’s orbit 

Sine of the eccentric anomaly of one of the transfer orbits 

Hyperbolic sine of the eccentric anomaly of one of the transfer orbits 

Number of Julian centuries since the year 2000 AD 

TIM minus half a second 

TIM plus half a second 

True anomaly of Earth or Moon orbit 

True anomaly of one of the transfer orbits (Gamacalc Subroutine) 



TIEM 

TIEM1 

TIEM2 

TIEMS 

TIEMT 

TIM 

TIMP 

TTMX 

TTMJ 

TRAJ 

TRAJE 

T1 

T2 

VCIRE 

VCIRX 

VEL 

VELE 

VELE2 

VLAT 

VLON 

VMAG 

VPE 


Iterated Earth (or Moon)-to-LP time of flight (seconds) 

Earth or Moon-to-LP time of flight guess (seconds), using VELE1 
Earth or Moon-to-LP time of flight guess (seconds), using VELE2 
Temporary storage for Earth-to-LP time of flight 
Total time of flight 

Time of arrival at LP from Earth (or Moon), in centuries since the year 2000 AD 
Time now 

Time of flight from Earth or Moon perigee to LP (Gamacalc Subroutine) 

Earth departure Julian date (where January 1, 2000 is day 2,451,545. Refer to 
Section C of The Astronomical Almanac of the Year 19881 

Text that describes trajectory as hyperbolic or elliptical (Gamacalc Subroutine) 

Text that describes trajectory as hyperbolic or elliptical 

One-half second before TIM 

One-half second after TIM 

Velocity of Earth circular orbit 

Velocity of the Earth or Lunar circular orbit (Gamacalc Subroutine) 

Velocity at LP of one of the transfer orbits (Veltrans Subroutine) 

Velocity at LP of the Earth (or Moon)-to-LP trajectory 
One foot per second more than VELE 
Latitude component of the velocity vector 
Longitude component of the velocity vector 
Velocity vector magnitude 

Perigee velocity of the departure body-to-LP trajectory 

Perigee velocity for one of the transfer orbits (Gamacalc Subroutine) 


VPX 



w 

VPMIN 

VXE 

VXLA 

VXLO 

VXR 

VXX 

VYE 

VYLA 

VYLO 

VYR 

VYX 

VZE 

VZLA 

VZLO 

VZR 

VZX 

XDLO 

XXBC 

XXL 

XXLP 

YDLO 

YYL 

YYLP 


Velocity at LP of one of the transfer orbits (Gamacalc Subroutine) 

Minimum velocity at LP for Earth-to-LP trajectory with maximum flight time 
X-coordinate of velocity vector at LP for Earth (or Moon)-to-LP trajectory 
X-component of the latitude component of the velocity vector 
X-component of the longitude component of the velocity vector 
X-component of the radial component of the velocity vector 
Total X-component of the velocity vector 

Y-coordinate of velocity vector at LP for Earth (or Moon)-to-LP trajectory 
Y-component of the latitude component of the velocity vector 
Y-component of the longitude component of the velocity vector 
Y-component of the radial component of the velocity vector 
Total Y-component of the velocity vector 

Z-coordinate of velocity vector at LP for Earth (or Moon)-to-LP trajectory 

Z-component of the latitude component of the velocity vector 

Z-component of the longitude component of the velocity vector 

Z-component of the radial component of the velocity vector 

Total Z-component of the velocity vector 

X-coordinate of the velocity of the Moon 

Distance from Earth geometric center to baricenter 

Array of distances in the X-direction from the Earth to each LP 

Distance in the X-direction from the departure body (Earth or Moon) to the LP’ 
X-coordinate 

Y-coordinate of the velocity of the Moon 

Array of distances in the Y-direction from the Earth to each LP 

Distance in the Y-direction from the departure body to the LP’s Y-coordinate 


zz 

ZZLP 


Distance in the Z-direction from the Earth-Moon plane to the LP (always 0) 
Distance in the Z-direction from the Earth-Moon plane to the LP (always 0) 



Appendix D. Detailed Program Description 


Librate Main Program 

1. Declare matrices (XXL, YYL, PAGE1, PAGE2, PAGE3, AINCEP, FLTIMP). 

2. Open the output file (LIBRATE.OUT). 

3. Define the program constants. 

4. Record the current date and time (DATP and TIMP). 

5. Read the program inputs: 

a. BODY (EARTH or MOON) 

b. HPE (perigee altitude of Earth orbit) 

c. NLP (libration point number) 

d. FLITM (flight time constraint) 

6. Calculate RPE, the distance from the Earth’s center to Earth perigee orbit. 

7. Calculate RM, Earth-Moon distance (Call Moon Subroutine). 

8. Calculate YDLO, the orbital velocity of the Moon (Call Posvelmo Subroutine). 

9. Calculate RREM, Earth-Moon distance in feet. 

10. Calculate libration point locations, maximum flight time (AMAXFT), minimum velocity 
at libration point (VPMIN) — In-program subroutine. If flight time input (FLT1M) is 
greater than the maximum flight time (AMAXFT) then re-input flight time (step 5). 

11. Update libration point locations, maximum flight time (AMAXFT), minimum velocity at 
libration point (VPMIN) — In-program subroutine. 

12. Save current time of flight estimate (T1EMS) 



13. Calculate time of flight from perigee to LP (T1EM1), cosine of the flight path angle 
(COSGAMX) at LP, velocity at perigee needed to reach LP (VPX), velocity for circular 
orbit at perigee (VCIRX) from current values of LP distance (RRL), and velocity at LP 
(VELE) — Call Gamacalc Subroutine. 

14. Increment velocity at LP (VELE) by 1 ft/sec to get VELE2. 

15. Calculate another time of flight (TIEM2), cosine of the flight path angle (COSGAMX), 
velocity at perigee needed to reach LP (VPX), velocity for circular orbit at perigee 
(VCIRX) from current values of LP distance (RRL), and velocity at LP (VELE2) — Call 
Gamacalc Subroutine. 

16. Estimate a delta-v needed to be added to VELE satisfy the time of flight constraint using 
a linear extrapolation of the two values determined in steps 13 and 15 (DEL VELE). 

17. If magnitude of DEL VELE increment is larger than 500 ft/sec than limit magnitude of 
DEL VELE to 500 ft/sec. Add the velocity increment to VELE. 

18. If the new velocity at the LP (VELE) is smaller than the previously calculated minimum 
velocity (VPMIN) then set VELE = VPMIN. 

19. Calculate new time of flight (TTEM), cosine of the flight path angle (COSGAMX), 
velocity at perigee needed to reach LP (VPX), velocity for circular orbit at perigee 
(VCIRX) from current values of LP distance (RRL), and velocity at LP (VELE) -- Call 
Gamacalc Subroutine. 

20. If the time of flight is zero then go to output section of program. 

21 . Total time of flight (T1EMT) equals time of flight from Earth or Moon (TIEM). 

22. If calculated flight time (ITMEE) is not within 1 hr. of the desired time then iterate the 
flight times again starting at step 11. 



23. Determine the flight path angle from the cosine of the angle (COSGAMX). 

24. Calculate longitude of libration point with respect to -X axis. 

25. Calculate latitude of libration point with repect to Earth-Moon line. 

26. Calculate azimuth angle of trajectory at LP. 

27. Convert velocity at LP to rectangular coordinates (Call Veltrans) 

28. Calculate the velocity increment needed to adjust the velocity heading at the LP such that 
the resulting orbit is stationary with respect to the Earth-Moon line (DELVEL) — angular 
velocity at LP, with respect to the Earth, is equal to that of the Moon (OMEGM). 

29. Determine velocity increment needed at perigee of initial departure orbit to attain an orbit 
containing the LP (DVCIRE). 

31. Save previous delta-v total as DVTSAV. 

32. Calculate total delta-v (DVTOTAL) from DELVEL + DVCIRE. 

33. Calculate position and velocity of the Moon: Determine exact Julian date and Call 
Posvelmo. 

34. Update distance from Earth to Moon in feet (RREM). 

35. If time of flight (HEM) is not within 0.09 hours of the time of flight before the last 
iteration (HEMS) then go back to step 11 for more iterations. 

36. Save current inclination, flight time, total delta-v, heading correction AV, and departure 
AV in the arrays AINCEP, FLTTMP, PAGE1, PAGE2, and PAGE3, respectively. 

37. Decrement flight time (FTIM) by 2 ft/sec. 

38. Go back to step 12 to calculate the trajectory for the new flight times unless the entire 
column has been completed (19 flight times for each inclination). 

39. Reset flight time back to input value (FLTIM). 



40. 


Increment inclination (AINCE) by 10°. 

41. Go back to step 12 to calculate the trajectory for the next column until the inclination has 
reached 90°. 

Produce three pages of output in the file LIBRATE.OUT using the saved arrays in step 
36 via Subroutine LPPAGE. 


42. 



Appendix E. Detailed In-Program Subroutine Description 

Librate In-Program Subroutine 

1. Update values for Earth-Moon position in feet (RREM), Earth-Moon baricenter (XXBC), 
and the angular velocity of the Moon (OMEGM). 

2. Update all libration point positions with respect to Earth-Moon line (XXL, YYL, ZZLP) 

3. Update libration point position in use (XXLP, YYLP), and correct values if departing 
from the Moon. 

4. Update distance from departure body to libration point (ERL), and longitude of LP with 
respect to departure body. 

5. Calculate minimum velocity (VPMIN) at libration point for orbit with apogee RRL and 
perigee RPE. 

6. If current iterated velocity at the LP (VELE) is less than minimum velocity then reset 
VELE to VPMIN. 

7. Call Gamacalc Subroutine to determine the maximum flight time (AMAXFT) for an 
orbit containing RRL and RPE, and having an apogee velocity of VPMIN. 

8. If input value for flight time (FLTTM) is greater than the maximum flight time then print 
the maximum flight time allowed (AMAXFT) and go back to step 5 to re-input the flight 
time constraint. 

9. Return 



Appendix F. Subroutine GAMACALC 


The subroutine GAMACALC receives the parameters orbital perigee radius (RPX), orbital 
velocity at LP (VV), orbital radius at LP (RRX), and the gravitational parameter of the body 
being orbited (CMUX). It calculates and returns time of flight (TIMX), the cosine of the flight 
path angle (COSGAMX), velocity at periapses (VPX), circular velocity at periapses (VCIRX), 
true anomaly (THETAX), and an indicator describing whether the orbit is elliptical or hyperbolic 
(TRAJ$). 

1. Initialize indicators to presume an elliptical orbit. (IHYPER, TRAJ$). 

2. Calculate the vis-viva parameter 

QQX = (RRX * VELX A 2) / CMUX. 

3. If the orbit is just barely hyperbolic (QQX is within one-millionth of 2), force the 
calculation to consider it elliptical (reduce QQX to just under 2). 

4. If the orbit is stil hyperbolic, reset the indicators to show this. Print a message on the 
screen announcing a hyperbolic orbit. 

5. Calculate the semi-major axis of the orbit 

AAX = RRX / (2-QQX). 

6. If the semi-major axis is very large (greater than 10 A 12) or very small (less than -10 A 12), 
the orbit is trapped near a parabolic trajectory. Make it hyperbolic: 

AAX = -10 A 12. 

7. Calculate the orbit eccentricity 

EEX = 1 - (RPX/AAX). 

8. Calculate the semi-latus rectum 

PPX = AAX * (1 - EEX). 

9. Calculate the flight path angle 

a. The angular momentum of the orbit at perigee is 

HP = RPX * VELPERIGEE * cos(GAMAX). 

But at perigee, GAMAX is zero, so 
HP = RPX * VELPERIGEE. 

b. At LP, angular momentum is 

HX = RRX * VELX * cos(GAMAX). 

c. Angular momentum is constant along a given orbit, so 

HP = HX 

RPX * VELPERIGEE = RRX * VELX * cos(GAMAX) 

GAMAX = arccos((RPX * VELPERIGEE) / (RRX * VELX)) 


d. The velocity at perigee is 

VELPERIGEE = sqr((CMUX * RAPOGEE) / (AAX * RPX)). 

Since RAPOGEE / AAX = (1 + EEX), then 

VELPERIGEE = sqr (CMUX * (1 + EEX) / RPX) or 
RPX * VELPERIGEE = sqr(RPX * (1 + EEX) * CMUX). 

e. Substituting (d) into (c) above yields 

GAMAX = arccos((sqr(RPX * (1 + EEX) * CMUX) / (RRX * VELX)) 

f. Substituting from (2) above yields 

GAMAX = arccos(sqr((RPX * (1 + EEX)) / (RRX * QQX))). 

10. Calculate the perigee velocity 

VPX = sqr(CMUX * (1 - EEX) / RPX). 

1 1 . Calculate the circular velocity at perigee 

VCIRX = sqr(CMUX / RPX). 

12. Calculate the true anomaly 

THETAX = arccos(((PPX / RRX) - 1) / EEX) 

13. Calculate the eccentric anomaly 

AEX = arccos((EEX + cos(THETAX)) / (1 + EEX * cos(THETAX))). 

14. Compare the orbital radius at LP (RRX) to the apogee radius (AAX * (1 + EEX)). If the 
orbital radius at LP is greater than the apogee radius, print a message on the screen 
indicating the difference in nautical miles. Also display the following: 

a. LP latitude (ALAT) 

b. LP longitude (ALON) 

c. vis-viva parameter (QQX) 

d. semi-major axis (AAX) 

e. eccentricity (EEX) 

f. flight path angle (GAMAX) 

g. perigee velocity (VPX). 

15. Calculate the time of flight. 

a. If the orbit is elliptical: 

T1MX = sqr(AAX A 3/CMUX) * (AEX - EEX * sin(AEX)). 

b. If the orbit is hyperbolic: 

TIMX = sqr(-AAX A 3/CMUX)*(EEX*sinh(EEX) - log(cosh(EEX) + 
sinh(EEX))). 


Appendix G. Subroutine POSVELMO 


The subroutine POSVELMO receives the parameter TIM (number of Julian centuries from the 
year 2000) and returns the moon’s position and velocity at that time. Specifically, it returns the 
moon’s distance from the Earth’s center, in Earth radii (RRM); velocity in the x-direction (along 
the Earth-Moon line), in feet per second (XDLO); and velocity in the y-direction (direction of 
the moon’s orbit), in feet per second (YDLO). 

1. Calculate a fraction of time that represents half a second. 

DELT = 0.5 seconds / (36525 days/century * 24 hrs/day * 3600 seconds/hr) 

= 1.58440E- 10 centuries. 

2. Call the subroutine MOON with the parameter (TIM minus DELTA) to determine the 
moon’s distance in Earth radii at half a second before TIM. This distance is RM1 . 

3. Call the subroutine MOON with the parameter (TIM plus DELTA) to determine the 
moon’s distance in Earth radii at half a second after TIM. This distance is RM2. 

4. Calculate the velocity of the moon in the -x (radial) direction. 

XDLO = [(RM2 Earth radii - RM1 Earth radii) / (1 sec)] * 20,925,741 ft/radii. 

5. Calculate the average radius of the Lunar orbit during the one second centered on TIM. 

RRM = (RM1 + RM2) / 2. 

6. Determine the radius of the Lunar orbit from the Earth-Moon barycenter. 

RRMB = RRM - 0.7412789 Earth radii. 

7. Calculate the moon’s velocity in the direction of its orbit. 

a. Moon’s apogee (Apo) = 62.83308 Earth radii. 

b. Moon’s perigee (Per) = 55.68264 Earth radii. 

c. Eccentricity (e) = (Apo - Per) / (Apo + Per) = 0.06033. 

d. Semi-latus rectum (p) = Apo(l -e A 2) = 62.60439 Earth radii 

= 1,310,038,967 feet. 

e. Earth’s gravitational parameter (mu) = 1.407646822E+16 ft A 3/sec A 2. 

f. Angular momentum (h) = sqr(mu * p) 

= 4.29427E+12 ft A 2/sec * (1 earth radii / 20,925,672.57 ft) 

= 205,215.4 ft*Earth radii/second. 

g. Y- velocity (YDLO) = h / RRMB = 205,215.4/RRMB. 



Appendix H. Subroutine MOON 


The subroutine MOON receives the parameter T (number of Julian centuries from the year 2000) 
and returns the approximate location of the moon in geocentric coordinates at that time. 
Specifically, it returns the right ascension of the moon (RAM), declination of the moon 
(DECM), and distance to the moon in Earth radii (RM). The formulae are from The Astronomi- 
cal Almanac of the Year 1984 . page D46. All degrees are converted to radians with the 
conversion factor C = tc/180. 

1 . Calculate the ecliptic coordinates of the moon. 

a. Moon’s ecliptic longitude 

LAM = 218°. 32 + 481267°.833T 

+ 6°.29 * sin(134°.9 + 477198°.85T) 

- 1°.27 * sin(259°.2 - 413335°.38T) 

+ 0°.66 * sin(235°.7 + 890534°.23T) 

+ 0°.21 * sin(269°.9 + 954397°.70T) 

- 0°.19 * sin(357°.5 + 35999°.05T) 

- 0°.U * sin(186°.6 + 966404°.05T) 

b. Moon’s ecliptic latitude 

BETA = 5°. 13 * sin( 93°.3 + 483202°.03T) 

+ 0°.28 * sin(228°.2 + 960400°.87T) 

- 0°.28 * sin(318°.3 + 6003°.18T) 

- 0°.17 * sin(217°.6 - 407332°.20T) 

c. Horizontal parallax 

PIE = 0°.9508 

+ 0°.0518 * cos(134°.9 + 477198°.85T) 

+ 0°.0095 * cos(259°.2 - 413335°.38T) 

+ 0°.0078 * cos(235°.7 + 890534°.23T) 

+ 0°.0028 * cos(269°.9 + 954397°.70T) 

d. Semi-diameter of moon’s orbit 

SD = 0.2725 * PIE 

e. Distance to the moon in Earth radii 

RM = 1 / sin(PIE) 


2. Form the geocentric direction cosines to rotate into geocentric coordinates. 

a. 1 = cos(B ET A)cos(LAM) 

b. m = 0.9175*cos(BETA)sin(LAM) - 0.3978*sin(BETA) 

c. n = 0.3978*cos(BETA)sin(LAM) + 0.9175*sin(BETA) 

where 1 = cos(DECM)cos(RAM), m = cos(DECM)sin(RAM), n = SIN(DECM). 



3. 


Then: 


a. RAM = arctan(m/l) [right ascension] 

b. DECM = arcsin(n) [declination] 

The errors will rarely exceed 0.2 Earth radii in distance (RM), 0.3° in right ascension (RAM), 
and 0.2° in declination. 



Appendix I. Subroutine VELTRANS 


The subroutine VELTRANS converts an orbital velocity vector into rectangular coordinates (see 
figure J-l). Parameters received by the subroutine are velocity (VEL), flight path angle 
(GAMA), azimuth (AZM), latitude above the Earth-Moon plane (ALAT), and longitude from 
the Earth-Moon line (ALON). A set of intermediate calculations is performed to express the 
velocity vector in terms of a radial component, a latitudinal component, and a longitudinal 
component. Each of these three components is further resolved into x-, y-, and z-components. 
Finally, all three x-components, all three y-components, and all three z-components are summed 
to provide the total x-, y-, and z-components of velocity. 

1 . Conversion of velocity vector into spherical coordinates. 

From the geometry, the radial component of velocity, R, is calculated to be 
VEL * sin(GAMA). (See/igure J-2). 

The component along the orbital path, Rif, is 
VEL * cos(GAMA). 

This orbital path component of velocity is further resolved into a latitude component, 

LAT, and a longitude component, LON (see figure J-3). Again, from the geometry, 

LON = Ri* sin(AZM) and 
LAT = Ri* cos(AZM). 

2. Conversion of spherical coordinates into rectangular coordinates. 

a. Conversion of radial component into rectangular coordinates. 

Refer to figure J-4. The projection of R onto the x-y plane is 
k * cos(LAT). 

The negative x-component of this is 
R * cos(LAT) * cos(LON) 
so the x-component, X, is 

-R * cos(LAT) * cos(LON). 

The negative y-component of R is 
k * cos(LAT)* sin(LON) 
so the y-component, Y, is 

-R * cos(LAT) * sin(LON). 

The z-component, 2, is 
k * sin(LAT). 

b. Conversion of latitude component into rectangular coordinates. 

Refer to figures J-5 and J-6. LAT is perpendicular to the radial vector, R. A line 
in the z-direction that meets the tip of LAT and intersects the radius vector 
produces the angles a and b, where 
b = 90 - LAT and 
a + b + 90 = 180. Therefore, 
a = LAT. 

From the geometry, the z-component of LAT, ZLAT, is 
LAT * cos(LAT). 



The projection of LAT onto the x-y plane is 
LAT * sin(LAT). (see, figures J-7). 

The x-component of this, XLAT, is 
LAT * sin(LAT) * cos(LON). 

The y-component of this, YLAT, is 
LAT * sin(LAT) * sin(LON). 

c. Conversion of longitude component into rectangular coordinates. 

« 

Refer to figures J-8 and J-9. LON is always parallel to the x-y plane, so the z- 
component of LON, ZLON, is always zero. Using the same proof as in (b) above, 
it can be seen that the angle between LON and the y-component of LON is equal 
to LON. From the geometry, the x-component of LON, XLON, is 
L6N * sin(LON). 

The negative y-component of LON is 
LON * cos(LON), 
so the y-component, YLON, is 
-L6 n * cos(LON). 

Sum of the rectangular coordinates. 

All of the x-, y-, and z-components are summed to provide the complete rectangular 

coordinates of the velocity vector. 

VXX = X + XLAT + XLON 
YYX = Y + YLAT + YI^ON 
VZX = Z + ZLAT + ZLON. 
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Figure J-6 
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